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We analyze actual methods that generate smooth frame fields both in 
2D and in 2>D. We formalize the 2D problem by representing frames as 
functions (as it was done in 3D), and show that the derived optimization 
problem is the one that previous work obtain via “representation vectors”. 
We show (in 2D) why this non linear optimization problem is easier to 
solve than directly minimizing the rotation angle of the field, and observe 
that the 2D algorithm is able to find good fields. 

Now, the 2D and the 3D optimization problems are derived from the 
same formulation (based on representing frames by functions). Their ener¬ 
gies share some similarities from an optimization point of view (smooth¬ 
ness, local minima, bounds of partial derivatives, etc.), so we applied the 
2D resolution mechanism to the 3D problem. Our evaluation of all exist¬ 
ing 3D methods suggests to initialize the field by this new algorithm, but 
possibly use another method for further smoothing. 

Categories and Subject Descriptors: 1.3.5 [Computational Geometry and 
Object Modeling]: Curve, surface, solid, and object representations 

General Terms: Frame field, 3D mesh 
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1. INTRODUCTION 

In computer graphics, a frame field can be defined on a surface 
(2D) or inside a volume (3D). For each point of the domain, it 
defines a set of 4 (resp. 6) unit vectors invariant by rotations of 7r/2 
around the surface normal (resp. around any of its member vector). 
The main motivation to study these fields is to split the quad and 
hexahedral remeshing problems into two steps: (1) the design of a 
smooth frame field, (2) and the partitioning of the domain by quads 
or hexes aligned with the frame field. Our objective is to unify the 
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formulation of the 2D and 3D frame field design problems, and to 
use it to extend an efficient 2D algorithm to the 3D case. 

In most cases, frame field design is formalized as the follow¬ 
ing optimization problem: find the smoothest frame field subject 
to some constraints. What makes them different from each others 
is obviously the dimension of the frames (2D or 3D), but also the 
definition of the field smoothness, the expression of the constraints, 
and the optimization method. Interestingly, the 2D case and the 3D 
case are addressed by very different strategies: 

—In 2D, the frame field design problem can be restated as a vector 
field design problem thanks to the introduction of the “represen¬ 
tation vector”. In local polar coordinates, each vector of a frame 
has the same angle modulo 7r/2, if we multiply it by 4 we obtain 
a unique representation vector (modulo 27r). It is easy to derive 
optimization algorithms acting on the representation vectors. For 
simplicity reasons, we limit ourselves to planar frame fields and 
use th e algorithm proposed by Kowalski et al |Kowalski et al.| 
|20I2| as reference. 

—In 3D, it is not possible to extend the idea of “r epresentation vec¬ 
tor”. Instead, Huang et al [Huang et al. 20II| propose to repre¬ 
sent frames by functions defined on the sphere, refer to figure 
for an illustration. A definition of the field smoothness is derived 
from this representation and optimized in a two step procedure: 
(1) initializati on based on optim ization of spherical harmonics 
coefficients in [Huang etal. 2011| or front propagation of bound¬ 
ary constraints in [Li et al. 2012| , followed by (2) smoothing it¬ 
erations performed by L-BFGS on Euler angles representation of 
frames. 

Thus our goal is to better understan d how 2D and 3D pr oblems 
are related to each o ther and to extend [Kowalski et al. 201^ to 3D. 
We first show that [Kowalski et al. 20I2| can be derived with the 
formalization approach inherited from the 3D case, and then we 
extend it to 3D by the same logical flow. Busy readers interested 
in on ly rep roduction of the method can skip to implementation sec¬ 
tion §5.5] the only required tool is a linear solver, all calculations 
are given explicitely. 

The 2D algorithm with frames represented by 
functions 

Solutions developed for 3D are very different from 2D solutions 
because the “representation vector” trick does not extend nicely 
into 3D. To unify both problems, we propose to go in the other di¬ 
rection we apply the functional frame representation to the 2D 
case. By doing so, we found another way to introduce the “repre¬ 
sentation vectors”: they appear as coefficie nt ve ctor of the function 
decomposed in the Fourier function basis §2.2| Following the logi¬ 
cal flow introduce d for the 3D case, we derive an estimation of the 
field smoothness §2.3| and formulate the corresponding optimiza- 
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Fig. 1. Orthonormal frames (close-up (a)) are represented by spherical harmonic functions (close-up (b)), attached to each vertex of a tetrahedral mesh. 
Streamlines and singularities of the field are shown in yellow and red, respectively. 


tion problem §2.6| We end up with exactly the formulation of 2D 
methods based on the “representation vectors”. 

This common formulation of the 2D problem is not the simplest 
one: the difference between adjacent frames is not evaluated by the 
rotation between them as in IBommes et al. 2013||Ray et al. 2008| , 
but approximated by the Euclidian distance in a plane. We obs erve 
the impact of this approximation on the objective funct ion §2.7| and 
show how it simplifies the optimization problem §2.8| 

The resolution mechanisms of the 2D optimization problem 
are strongly inspired by the geometric intuition of the represen¬ 
tation vector field. To make abstraction of this int uition, we re¬ 
explain the algorithm from [Kowalski et al. 2012| with the no¬ 
tations/vocabulary introduced by the functional representation of 
frames. 


Extension of the 2D algorithm to 3D 

Now that we can describe the 2D algorithm without referring to 
the representation vector, it is possible to extend it to 3D. Using the 
notations introduced in the 2D case, we describe the 3D version in 
^and extend the 2D optimization mechanism to 3D §3.5| 

A first difficulty was to find the expression of the boundary 
conditions because the b oundary frames are free to rotate around 
the surface n ormal §3.4| Incompl ete enforcing of this condition by 
Huang et al [Huang et al. 201 1| results in a poo r initialization of 
the optimization procedure as evaluated in §3.6.1| 

The second difficulty is that the frame is represented in a 9D 
function basis, but the set of functions that corresponds to a frame 
has dimension 3 (the 3D rotation that brings the axis aligned frame 
to the current frame). T he extension of the no rmalization of the 
representation vector in [Kowalski et al. 20T^ becomes: find the 
nearest point on the 3D manifold of the set of functions represent¬ 
ing a frame. _ 

Our extension of [Kowalski et al. 201^ nice ly completes the tool 
set of 3D frame field design algorithms §3.6| Our initial field has 
lower energy, often a better topology, and is more robust to sur¬ 
face noise. The optimization step is easy to implement from the 
initialization step, is competitive when the initial topology is good 
enough, but does not outperform the current state of the art. 


Previous works 

The orientation of objects is commonly represented by symmetric 
tensors in physics to model the orientation distribution of fibers. For 
example, Moakher [Moakher 2009] introduced the notion of cubic 
orientation distribution functions. However, in computer graphics, 
more compact representations are often preferred: “representation 
vectors” in 2D and vectors of spherical harmonics coefficients in 
3D. 

On surface. The optimization of a frame field inside a 2D shape 
is very similar to the optimization of direction fields on surfaces [j 

The problem of direction field design on surfaces was introduced 
by [Hertzmann and Zorin 2000] for orienting strokes in non photo 
realistic rendering. Directions are represented by an angle rota¬ 
tion 0 per vertex, and the smoothing is performed by a non linear 
solver (BEGS). The ’’representation vector” was not introduced yet, 
and the optimization mechanism was similar to actual approach of 
3D frame fields smoothing. Solving with a repr esentation vector 
V = {r cos{0), r sm{0)) for each 0 was used in [Ray et al. 2006) 
for faster results, and impr oved later for better control over the field 
topol ogy [Ray et al. 2009|. Based on the similar representative vec¬ 
tors, [Palacios and Zhang 2007a| propose to control the field topol¬ 
ogy by local operations. For direction fields without constraints, 
[Knoppel et al. 201^ proved that the norm of the representative 
vector does not affects the result, leading to finding optimal direc¬ 
tion fields by solving an eigenvector problem. 

Directly working with angle 6 allows to perfectly control the 
field topology, but at the expense of solvi ng a mixed integer sys¬ 
tem [Ray et al. 2008l|Bommes et al. 200^ . In [Palacios and Zhang] 
|2007a| , representative vectors are introduced with its duality with 
N'^n order traceless symmetric tensors. This relation is very useful 
to unify 2D and 3D frame fields. 

Inside a volume. The pioneer work of [Huang et al. 2011| dis¬ 
covered that each frame can be represented by a spherical har- 


^The differences between these two problems (angle defect, hard con¬ 
straints, curvature fitting term, etc.) have an impact on the optimization 
algorithm, as detailed in the supplemental material 
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monic. Their initialization step defines a smooth spherical har¬ 
monic field then, for each sample, defines the frame that better 
aligns with the spherical harmonic. This initial solution is then im¬ 
proved by rotating each frame. These rotations are defined by Eu¬ 
ler angles and obtained by minimizing the field smoothness with 
L-BFGS a nd enfo rcing the boundary alignment by a penalty term. 

Li et al |2012| propose an alternative initialization method. They 
optimize ?l 2D frame field on the volume boundary, convert it to a 
3D frame field by adding the surface normal and its opposite, then 
propagate it inside the volume. The resulting fi eld is then smoothed 
by optimizing a rotation for each frame, as in |Huang et al. 2011) , 
but with an improved optimization scheme. They also optimize the 
singulari ty graph of the fie ld by local combinatorial operations, as 

done by [Jiang et al. 2014|. _ 

Our extension of [Kowalski et al. 2012| is optimizing for the 
same objective function, but with a very different solution mech a- 
nism. It performances are compared with previous work in §3.6| 

2. FUNCTIONAL REPRESENTATION OF FRAMES 
IN 2D 

This section introduces how to optimize 2D frame fields using a 
functional representation for each frame. While we do not claim 
any technical contribution in this section, we think that it is impor¬ 
tant to reformulate existing concepts using the functional represen¬ 
tation, because 3D case inherits exactly the same difficulties and 
the intuition we gained in 2D helps to motivate the choices made 
for 3D fields. We derive an energy and the boundary conditions 
from this new representation. The resulting optimization problem is 
exactly the one usually solved by direction field algorithms based 
on the standard “representation vector”. We show on a toy and a 
real examples that this optimization problem is much easier than 
directly minimizing the frame rotation. We then present the algo¬ 
rithm [Kowalski et al. 201^ that we extend to 3D in the following 
section. 

2.1 Problem settings 

Given a 2D shape, frame field design in 2D consists of finding a 
smooth frame field aligned with the boundary of the shape. We for¬ 
mulate it as minimizing the field curvature, based on the following 
definitions 0 

—A frame is a set of 4 unit vectors / = {fi},i G [0,3] that is 
invariant by a rotation of tv/2 (Figured. It can be represented 
by the angled such that Vi,/i = (cos(^i7r/2), sin(^+i7r/2)). 
—A frame field is a frame per vertex of a 2D shape triangulation. 
—The boundary constraint: a frame located on a boundary vertex 
must have one of its member vectors equal to the normal on the 
boundary. 

—The rotation angle between two frames is the angle AO of the 
rotation that transforms one frame into the other. This angle be¬ 
ing defined modulo 7r/2, we choose the AO with minimum ab¬ 
solute value. 

—The curvature of a frame field is the sum over each edge of 
the squared rotation angle between frames that are defined at the 
edge extremities. 

—A triangle is said to be singular0if the sum of the rotation angles 
over his boundary is not equal to 0. 


^The problem is directly presented in discrete settings. Interesting results on 
its (non trivial) continuous counterpart are given in supplemental materials. 
^Frame field topology is further discussed in supplemental material. 



Fig. 2. A 2D frame is a set / of 4 vectors /o, /i, / 2 , /s invariant under 
rotation by tt/ 2. Its angle representation is the rotation 0 between the global 
axis X and /o. 

Representing frames by angles is simple, but it makes the field 
curvature hard to optimize [Bommes et al. 200^|Ray et al. 2008) , 
and it does not extend nicely in 3D. For these reasons, we pro¬ 
pose an alternative representation based on functions, and use this 
curvature based formalization as a reference. 


2.2 Functional approach: frame representation 

The frame aligned with coordinate axes is called the reference 
frame / = {(1,0), (0,1), (—1,0), (0, —1)}. Instead of using the 
angle approach, we represent it by the function F{a) = cos(4a) 
with a G [0, 27r[ (Figure[^left), that exhibits the same 7r/2 rotation 
invariance as the frame. 

Any other frame / can be obtained by a rotation of / by angle 
0. The functional counterpart is to rotate the graph of the function 
F, namely any frame / can be represented by a function F(a) = 
F{a — 0) — cos(4(a — 0)) with a G [0, 27r[ (Figure [fright). 

A compact representation of these functions is given by the 
trigonometric relation/(a) = cos(4a —4^) = cos(4^) cos(4a) + 
sin(4^) sin(4a): we see that a frame function F can be rep¬ 
resented by a 2D vector of coefficients a — (ao,ai)^ = 
(cos(4^), sin(4^))^ in Fourier basis B = (cos(4a), sin(4a)) i.e. 
F^Ba. 

A coefficient vector a is feasible if and only if there exists 0 
such that a — (cos(4^), sin(4^))^. Geometrically, a is constrained 
to live on a curve parameterized by 0. This curve represents, in 
coefficient space, all possible rotations of the reference frame. In 
2D, this curve is the unit circle, so the constraint on a is simply : 
a^a = 1 . 

At this point, we can observe that the coefficient vector a is ex¬ 
actly the representative vector used in the direction field literature. 
We can also notice that expressed in Cartesian coordinates, our ref¬ 
erence frame function F is the polynomial — 3 restricted 

to the unit circle, thus it is also eq uivalent to the traceless symm etric 
4^^ order tensors manipulated in (Palacios and Zhang 2007b| . 


2.3 Functional approach: objective function 

We have defined the field curvature as the sum over each edge of the 
squared difference between frames at the edges extremities. In our 
formalization, the difference between two frames (at vertices i and 
j) is no longer the rotation angle, but the norm of the difference 
between the corresponding functions : fQ^(F^ (a) — F^(a))^da. 
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0 



F{a) = 
F(a - e) 



Fig. 3. Parametric plot of the reference frame representation F{a) (left) 
and an arbitrary frame F{a) (right). The plot is made using x{a) = 
(1 + F(a)) cos(q:) and y{a) = (1 + F(a)) sin(Q:) for a G [0, 27r[. It 
is easy to see that corresponding frames are aligned with maxima of the 
representation functions. 


It leads to the new objective function: 




(Fya) - F\a)fda 


(Ba^ - Ba^fda 


= B^Bda^ (a^ - a*) 

As the function basis B is orthogonal, and all functions are of norm 
so the expression simplifies to: 


E = TrY,\ 


( 1 ) 


2.4 Functional approach: constraints 

As discussed in §2.2| the first constraint is clearly thal ihe variables 
a* must be feasible (i.e. there exists a frame represented by a*). 

The second constraint is that frames of boundary vertices i must 
have one member vector equals to the normal direction. If 0^ de¬ 
notes the normal direction, the frame can be directly fixed by satis¬ 
fying two equations: 

4 = cos(46>^) (2) 

a\ = sin(4^*) 

However, as we already have the feasibility constraint = 

1, enforcing only one equation has the same effect: 

Qq cos(4^*) + a\ sin(4^*) = 1. (3) 


2.5 Toy example 

It is natural to ask the question: “Does minimizing our energy min¬ 
imizes the field curvature as well?’' 

Two frames /* and are represented by a* and , both located 
on the unit circle. The field curvature measures the circle arc length 
between them, whereas our norm is the chord length between 
a* and . 

From a practical point of view, we want to produce smooth fields, 
so most edges will have low curvature. In this case the objective 
function E is almost proportional to the field curvature. If, how¬ 
ever, two adjacent frames are not similar (e.g. they are close to sin¬ 
gularities), then the function E is smoother than the field curvature, 

ACM Transactions on Graphics, Vol. 28, No. 4, Article 106, Publication date: August 2009. 



Fig. 4. Top row: interpolation problem with two locked extremity frames. 
Bottom row, left: field curvature plot. Bottom row, right: our objective func¬ 
tion E. The plots are made as functions of the rotations 6 i , 62 of interpo¬ 
lated frames. Both functions share same local minima Pq and Pi . 



Fig. 5. Two minima for the toy problem shown in Figure]^ Pq turns the 
frames counterclockwise while Pi turns clockwise. Pq minimizes energy 
E and has better field curvature. 


making it easier for the optimization algorithm to move singulari¬ 
ties. 

Let us illustrate our intuition on a very simple interpolation ex¬ 
ample: a chain of four vertices having its extremities locked. The 
toy problem is therefore to find two frames interpolating the ex¬ 
tremity frames. 

If we represent two free frames by their angle 0i and O 2 , we 
can plot and compare the field curvature versus our objective func¬ 
tion E (Figure]^. The field curvature is not smooth (it is piecewise 
quadratic) and we can observe that there are two local minima. Our 
objective function is smooth, and has exactly the same minima on 
this example. Note that it could also have a smaller number of min¬ 
ima e.g. if the constrained frames are more similar. 

Figure]^ shows the frames corresponding to minima Pq and Pi: 
they differ by the sense of rotation. The point Pq minimizes objec¬ 
tive function E and has better field curvature. In next two sections 
we expose current state of the art approach to minimization of the 
objective function. 
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2.6 Implementation 

We have to minimize our objective function E (eq. |T}) with linear 
equality constraints on boundary vertices (eq. ^ or eq. ([^) and 
quadratic e quality constraints = 1 for each vertex i. 

We use IKowalski et al. 201^ ’s algorithm to solve this prob¬ 
lem. It finds an initial solution by relaxing the quadratic feasibility 
constraints on a* and finding the nearest feasible solution. Then it 
performs a number of smoothing iterations to ameliorate the solu¬ 
tion. In order to respect the feasibility, the quadratic constraints are 
linearized at each smoothing step. 

Initializcition. Here, we relax the feasibility constraints, so we 
need to minimize the quadratic function E subject to linear bound¬ 
ary constraints. To do it, we simply replace the linear constraints 
by a strong penalty term in the objective function, leading to a new 
quadratic function to optimize. This penalty method is very simple 
and sufficient in practice. 

More precisely, the new quadratic function is expressed in the 
form II AX—where A is a matrix, X our variable vector (X 2 i = 
ttQ and X 2 i+i = nj) and 6 is a vector. The system of equations 
AX — 6 = 0 is constructed line-by-line: 

—initial objective function E: for each edge ij, we add two equa¬ 
tions (eq. ([T])): 

V^{X2i - X2j) = 0 
aA'(-^2z+1 — X2j + l) — 0 

—boundary constraints: for each boundary vertex i, we add two 
equations (eq. 

\X 2 i = A cos 

\X2i+i = Asin4^\ 

where we set A = 100 in our experiments. 

From A and b, we just need to solve the linear system A^AX = 
A^b to obtain a minimizer of ||AX — 6|p. Then from X we can 
obtain an initialization of a* by projecting corresponding vectors 
on the set of feasible coefficients: 

a* ^ {X2i,X2i+lV/\\{X2i,X2i+l)\\. 

Smoothing iterations. Each smoothing iteration is similar to 
the initialization problem, except that we add to the objective func¬ 
tion a new quadratic penalty term that correspond s to a linear ap¬ 
proxim ation of the feasibility constraint as done in IKowalski et al.| 
|2012[ p. 6]. As before it is expressed by a new set of linear equa¬ 
tions when constructing A and b\ for each vertex z, we add one 
equation A(X 2 iaQ + X 2 i+ia\ — 1) = 0, where a* denotes the 

solution obtained at the previous iteration. _ 

To solve linear systems we use OpenNL library [Levy ) : it au¬ 
tomatically constructs A^ AX = A^b from the set of linear equa¬ 
tions and then solves it by the conjugate gradient method. 

2.7 Toy problem revisited 

This section explains o ur op timization approach on the toy problem 
already presented in § |2.5| As we mentioned before, at the initial¬ 
ization step we relax the constraints of feasibility of Unfortu¬ 
nately we can not plot the corresponding energy since without the 
constraints it becomes four-dimensional. 

Top row of Figurej^shows the solution of the initialization stage. 
Intuitively, we allow the points a* not to be on the unit circle. Hence 
and lie on the chord between locked points and a^. Second 



Fig. 6. First row: initialization stage solution. Second row: correspond¬ 
ing functional interpretation. Third and fourth rows: frames obtained by the 
projection of initialization a* and after smoothing iterations. 

row of Figure shows the corresponding functions and the third 
row gives the frames obtained by projecting points a* on the circle 
of constraints. 

Note that the initialization stage produces the correct sense of ro¬ 
tation (Figurej^. However it does not directly produce the optimum 
point Pq. The reason being that the initialization stage produces 
points and (before normalization) in the way that all three 
chord segments are equal: ||ao — ai|| = ||ai — <^ 2 || = \\ci 2 — cisW. 
But after projecting the points onto the feasible circle (3rd row of 
Figure the corresponding arc lengths are not equal. Therefore, 
we need a few smoothing iterations to reach the optimum (Fig¬ 
ure [5]— bottom row). 

2.8 Results 

Figure. 13-left shows the optimization of the field curvature by a 
gradient descent algorithm, initialized by an axis aligned frame 
field. We obtain a frame field that bends to match the boundary 
constraints, but its singularities remain on the boundary. The sys¬ 
tem is only able to reach the local minima corresponding to the 
initial field topology. The field curvaturej^is 34.21. 

Figure [^middle shows the optimization using only the smooth¬ 
ing iterations, initialized by an axis aligned frame field. We observe 
that smoothing iterations were able to slightly move the singular¬ 
ities away from the border and even to merge two singularities. It 
results in a better field curvature (equal to 31.41, 24.01 and 23.95 
after 10^, 10^ and 10® iterations). 

Figure [fright shows that the initialization step alone finds a so¬ 
lution with a simple topology and a lower field curvature (20.84). 
Smoothing iterations further decrease the field curvature (to 17.91). 

These observations suggest that our initialization step is manda¬ 
tory, and smoothing iterations further improve the result. However, 
it is important to notice that each iteration takes approximately the 
same time as the initialization step. 

Boundary constraints. When working only with feasible so¬ 
lutions a single equation (equation is sufficient to enforce each 


^The value of the field curvature is relevant only for comparing different 
fields on the same mesh. 


ACM Transactions on Graphics, Vol. 28, No. 4, Article 106, Publication date: August 2009. 















6 


N. Ray and D. Sokolov 



Fig. 7. Evaluation of 2D frame field optimization algorithms: singular triangles are highlighted in red. Compared algorithms are: (left) steepest descent of 
the field curvature, initialized with an axis aligned frame field, (middle) smoothing iterations with our objective function E initialized with axis aligned frame 
field after 10^, 10^ and 10® iterations (from left to right), (right) initialization step alone (left) and the initialization step followed by 10® smoothing iterations 
(right). 


boundary constraint. Using it for the initialization step is wrong: for 
example, a normal constraint of angle ^ = 0 forces — 1 but let 
tti free. As illustrated in Figure it can produce very bad frames 

fields. Therefore we must use two separate equations pjf _ 

The re exists a similar issue in 317: Huang et al |Huang et al.| 
|201H use a 317 boundary condition that is not sufficient for the ini¬ 
tialization step. It leads to a poor initialization for their smoothing 
algorithm, making it very slow, and getting possibly loc ked w ith a 
bad initial topology. This issue is discussed in details in §3.6| 



Fig. 8. Boundary constraints for global optimization: if we only use one 
constraint (eq|3> the initialization step finds a constant frame field on a 
parallelogram (left). It is therefore mandatory to lock the frames to get the 
proper boundary constraints (right). 


3. OPTIMIZATION OF 317 FRAME FIELDS 

Our objective is to extend to 377 the optimization algorithm pre¬ 
sented in previous section. 

In 277, our framework allows to retrieve the representation vec¬ 
tor that was the key to efficient optimization of direction fields. 

We now extend our framework to 377, find a generalization of 
this representation vector, express the boundary alignment condi¬ 
tion with respect to this representation, and derive the optimization 
algorithm. 

3.1 Problem settings 

The problem is to define, inside a 377 shape, a smooth frame field 
that is aligned with the boundary of the shape. We are working in 
discrete settings on a tet mesh. The problem to minimize the field 
curvature is defined as follows: 

—The reference frame / is the set of 6 unit vectors forming nor¬ 
mals of a cube aligned with coordinate axes (Figure]^. 

—A frame is the reference frame rotated by a 3 x 3 matrix R\ 
f = Rf. Note that multiplying a matrix by a set is a slight abuse 
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of notation, it means that we obtain a new set where each vector 
is rotated by the given matrix. 

—A frame field is the definition of a frame for each vertex of the 
tet mesh. 

—The boundary constraint: The frame of a boundary vertex must 
have one of its member vectors equal to the normal of the bound¬ 
ary. 

—The rotation angle between two frames is the minimal angle of 
rotation that brings one frame to the other. 

—The curvature of a frame field is the sum, over each edge, of 
the squared rotation angle between adjacent frames. 

—A tet is called singular if any of its triangles is singularjj The 
triangle ijk is singular if and only if R^^ x R^^ x 77^* ^ Id, 
where R^^ denotes the rotation matrix that brings the frame /* 
to the frame /7. 

3.2 Frame representation 

The reference frame / is represented by the function F = 
where Yi^rn is the real valued spherical har¬ 
monic of degr ee I and order m. These harmonics are sometimes 
called tesseral [Ferrers 1877| p. 74]. The list of harmonics of d e- 
gree 4 can be found in IGorller-Walrand and Binnemans 1996| p. 
239]). Function F is defined as R® ^ R, but we are interested by 
its restriction to the unit sphere ^ R. 

Any other frame / can be obtained as a rotation of / by a matrix 
77. It is represented by the function F(P) = F{RP), where P is a 
point of the unit sphere (Figure]^. 

Yi^rn forms a functional basis over the unit sphere with an in¬ 
teresting property: applying a rotation to a spherical harmonic of 
degree I produces another harmonic of degree /.Asa consequence, 
since we represent the reference frame by a sum of two spherical 
harmonics of degree 4, each frame function F can be represented 
in the basis B = (U 4 _ 4 , 14 ,- 3 , • • • ? 14 , 4 )- Using it, we can rewrite 
the expression for the reference frame function as P = Bd with 


®We can define what a singular tet is, but we are not able to characterize 
them by an equivalent of the index in 2D. This fact is discussed in the 
supplemental material. 
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F{RP) = BRea 

Fig. 9. A frame / is the reference frame / rotated by a 3 x 3 matrix R. 
The plots of the corresponding functions F and F are also rotated by R, 
and their coefficients vectors verify a = where is a 9 x 9 rotation 
matrix given in §A.1| 


a — (^0,0,0,0, 0,0,0, frame / = Rf 

can be represented by F = Ba, where a = Rsa with Rb b eing 
a 9 X 9 rotation matrix acting on coefficients space. Appendix ]A. 1| 
describes the construction of the rotation matrices. 

A feasible coefficient vector a is a vector that can be written 
as a = Rsa where Rb is a 9D rotation matrix that can be de¬ 
rived from a 3D rotation. Geometrically, a is constrained to be on 
a manifold of dimension 3 embedded in the 9D coefficient space. 

At this point we can consider the coefficient vector a as an ex¬ 
tension of the representative vector used in the direction field litera - 
ture. It is also the representation introduced in [Huang et al. 2011) . 

3.3 Objective function 

As in 2D, the objective function is the sum, over each edge ij, 
of the squared difference between frames located at the edges ex¬ 
tremities. In our formalism, the difference between two frames is 
not the rotation angle, but the norm of the difference between 
the corresponding functions: {F^ (a) — F'^ (a))^da. It gives the 

energy: 

n27r 

E = y2 - F\a)fda 

ij F 

Here, the function basis B is orthonormal, so the expression sim¬ 
plifies to: 

£ = ^||a^'-a'f (4) 

ij 


3.4 Constraints 

There are two types of constraints: each coefficient vector a* must 
be feasible, and boundary frames must have one vector aligned with 
the normal of the volume boundary. The first constraint is presented 
in the frame representation section, and will be enforced by our op¬ 
timizer in a dedicated “projection” step (the 3D counterpart of the 
normalization of a in the 2D case). Here we focus on the boundary 
constraint. 

Smooth vertex. We assume first that there is only one normal as¬ 
sociated with the vertex, it can be computed as the average normal 
vector of incident boundary triangles. 

Case 1: the normal is equal to the z axis. 



Fig. 10. A 3D frame field produced on a {2D) disk (Left) produces the 
2D frame field we could obtain from the 2D algorithm (Right). 

Let us first consider the case where the fixed vector (the surface 
normal) is the axis z. If we rotate F around z by angle 0, we obtain 
a = RbO. with Rb being a rotation around z. The simple structure 
of Rb together with the null coefficients of a gives the equation: 


a= I A/^sin46»,0,0,0, y^>0>0,0, y^cos46» 


(5) 


As done in the construction of coefficient vectors in the 2D 
case, we can get rid of the angle 0 by replacing it by a vector 


- = (co,ci) = {\f^ cos46», 


^ sin AO 


)■ 


a = W^( 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 , 0 )T 

+ co(0,0,0,0,0,0,0,0,1)^ 

+ Ci(l, 0,0, 0,0,0,0,0,0)^ 


( 6 ) 

(7) 

(8) 


With this equation, all frames having a vector equal to z can be 
represented by the 2D vector c. As in the 2D case it comes with a 
norm constraint: Cq + cf = ^. 

The variable c defines the rotation of the frame around the sur¬ 
face normal i.e. a 2D frame field. The optimization of this 2D 
frame field using c as variables is exactly what we did in 2D by 
introducing the coefficient vector a. Our 3D solution restricted to 
the object boundary is therefore almost [^equivalent to our 2D so¬ 
lution (Figure p^. 

Case 2: the normal is not equal to the z axis. 

To handle this more general case, we rotate the constraint. If we 
want the vector n to be preserved, we first compute a rotation that 
brings z axis to n. From this rotation, we compute the correspond¬ 
ing 9D rotation matrix Rb, and derive the constraints: 

(9) 

( 10 ) 
( 11 ) 

This expression of the normal constraint gives us a set of 9 linear 
equations per boundary vertex. It introduces two new variables Cq 
and Cl, and a quadratic constraint Dc— 5/12. 

Note As in the 2D case, the boundary constraint has a simpler 
expression [Huang et al. 2011) : a^i?B(0,0,0,0,1,0,0,0,0)^ = 

that is valid only if all a* are feasible. Consequently, it cannot 
be used safely during the initialization step (see Figure p^. 


a = W 0,0,0,1,0,0,0,0)^ 

+ coi^s(0,0,0,0,0,0,0,0,1)^ 

+ cii?B(l,0,0,0,0,0,0,0,0)^ 


®The boundary has curvature that was not assumed in our 2D frame fields. 
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Non smooth vertices. Frames of vertices located on hard edges 
have to conform to more than one normal. These vertices have mul¬ 
tiple normal constraints, we pick two normals that are almost or¬ 
thogonal, perturb them (by rotations around their cross product vec¬ 
tor) to make them orthogonal, and compute the rotation that brings 
X to the first normal, and y to the second normal. We compute the 
corresponding coefficient space rotation Rb and fix the frame co¬ 
efficient vector a* to Rsa. 

3.5 Implementation 

We have to minimize our objective function (eq. (|^) with linear 
equality constraints on boundary vertices eq. i9}, quadratic equality 
constraints c* • = 1 on boundary vertices, and the constraint that 

each a* is feasible. _ 

As in the 2D case (in §2.6| ), our minimization algorithm 
(Algo.[B is formulated as a series of least squares problems (mini¬ 
mize \^X — ^IP), where A and b are constructed without the fea¬ 
sibility constraint of a* at the first iteration (initialization), and with 
a linear approximation of it in the subsequent iterations (smoothing 
iterations). 

—Initialization: Our variable vector X must represent the repre¬ 
sentation vectors a* but also the c* variables introduced to ex¬ 
press the boundary constraint. To do so, we first reorder vertices 
to have boundary vertices first We can then organize the vari¬ 
able vector Why blocks: W[9z+c?] = a^, and W[9n^; + 2i+(i] = 
where riy is the number of vertices. 

As in the 2D case, the matrix A and the vector b are constructed 
iteratively by adding new equations for the objective function 
(algorithmic and the boundary constraints (algorithm [C- In our 
approach, we do not explicitly enforce the feasibility of c* (c* • 
= 5/7), but it will be indirectly respected by the feasibility 
of a\ 

The projection of a* on the set of feasible coefficient vectors is 
no longer a simple normalization. Instead we perform, for each 
vertex, a gradient descent (algorithm [C initialized by a. More 
precisely, starting with a we rotate our current frame gradually in 
order to minimize the distance between the current frame func¬ 
tion and the function to be projected. The gradient is evaluated 
by calculating the variation of the norm induced by infinites¬ 
imal rotation matrices with Euler’s angles. 

—Smoothing iteration: For the linearized feasibility constraint of 
a% we must also add 3 extra variables per vertex to our system. 
These variables account for the position in a local basis of the 
tangent space of the 3T> manifold of feasible The introduction 

of these constraints in the matrix A is detailed in algorithm]^ 

This frame field design algorithm can be implemented without 
being expert in spherical harmonics. We give explicit construction 
of matrices Rb , in the appendix. The system to solve 

A^AX = A^b is simply a linear system with a positive definite 
matrix. We use the OpenNL library | |Levy | because the system can 
be directly constructed from the equations (lines of A and elements 
of b). 

In order to keep the algorithm easy to read, we did not detail how 
to lock frames for vertices with multiple normal constraints. 


^It is possible to increase the performances by 30% by doing a Hilbert 
sort in the boundary vertices block, and another one in the free vertices 
block 
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Algorithm 1: Frame field optimization 

Input: 

—A tetrahedral mesh M with: 

—Uy vertices including ni vertices with normal constraint 
—edges 8 

—number of smoothing iterations N 
Output: A frame /* for each vertex i 


1 sort(Ad, 8); I I vertex i is a boundary vertex i < ni 


2 foreach I e 0 ... N do 

3 HA and b will be constructed iteratively 

4 create matrix A with 0 rows and Oriy + 2ni + Sriy 
columns; 

5 create vector b of size 0; 

6 add_smoothing_terms(Ad, A, b); 

7 add_normal_constraints(Ad, A, b); 

8 // add constraints only if we are in a smoothing iteration 

9 if / > 0 then 

10 I addJocal_optim_constraints({a*}, M, A, b); 

11 end 


12 H solve A^AX — A^b 

13 X ^ callJeast_squares_solver(A, b)\ 


14 

15 

16 

17 

18 

19 end 


H find the frame for each vertex 
foreach z < do 

^ X[9i...9iES]; 

(/% a*) ^ closest_frame(aA; 


Algorithm 2: add normaLconstraints 


Input: A tetrahedral mesh Ad, matrix A, row b 
Output: Modified matrix A and vector b 

1 // enforcing normal constraints by quadratic penalty 

2 foreach i <ni do 

3 estimate normal n at vertex z; 

4 find Euler angles (a, /3, 7 ) to rotate z-axis to rzj _ 

5 find 9 X 9 rotation matrix Rb ; // see appendixiAl 

6 Hq i — Rb X (1,0,0,0,0,0,0,0,0)^; 

7 H^^Rb X (0,0,0,0,1, 0, 0,0,0)T; 

8 hg i — Rb X (0,0,0,0,0,0,0,0,1)^; 

9 A ^ 100; // quadratic penalty multiplier 

10 foreach d G 0 ... 8 do 

11 create vector row\ 

12 rozu[9z + d] ^ A; 

13 rozu [9n^; + 2z + 0] 

14 row^^riy -f 2z + 1] 

15 A.add_row(rozz;); 

6.push( A ^ 7 / 12/14 [d] ) 


16 

17 

18 end 


Xho[d]- 
■ A/zsidj; 


end 


3.6 Results 

It is impossible to compare frame field design algorithms only from 
the images and results present ed in the state of th e art papers. First 
of all, our implementation of [Huang et al. 2011| has significantly 
better performances compared to what was presented in the orig- 
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Algorithm 3: add_smoothing_terms 

Input: A tetrahedral mesh Ad, matrix A, row b 
Output: Modified matrix A and vector b 

1 foreach ij G do 

2 foreach d G 0 ... 8 do 

3 create vector row; 

4 row^i + d] ^ 1; 

5 row[^j-\-d]i -1; 

6 A.add_row(roi(;); 

7 5.push(0); 

8 end 

9 end 


Algorithm 4: add locaLoptim constraints 


2 

3 

4 

5 

6 

7 

8 
9 

10 
11 
12 

13 

14 

15 end 


Input: Previous solution {a*}^<n^;» tetrahedral mesh Ad, 
matrix A, row b 

Output: Modified matrix A and vector b 
1 foreach i < riy do 

Cx ^ X a*; // appendix A.2 

Cj/ — E^ X 
Cz ^ E^ X a*; 

A ^ 100; // quadratic penalty multiplier 
foreach d G 0 ... 8 do 

create vector row; 
row[^i + d] ^ A; 
row^rix + 2n; + 3z + 0] 
roiu[9n^ + 2nz + 3z + 1] 
row^rix + 2n; + 3z + 2] 

A.add_row(rou;); 

6.push(Aa*[d]); 


-Acj;[d]; 

-Acy[d]; 

-Ac^[d]; 


end 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 


Algorithm 5: closestArame 
Input: 9-component vector q 
Output: A frame / and its representation vector a 

CL i — a; 

s ^ 10“^; // optimization step size 
e ^ 10“^; // step threshold 
q^q/\q\; 

while True do 

g ^ {q^E^a^ q^E^a, q^E^a); H gradient in point a 
if 11^II < £ then 
I break; 

end 

Rb ^ R%{s ■ g[0]) X Rl{s ■ g[l]) x • g[2])- 
R ^ R-{s ■ s[0]) X Ry{s ■ s[l]) X R%s ■ g[2]y, 

CL i — R^CL ; 

end 

return /, a; 


inal paper. Next, Li et al [Li et al. 2012| did not present any 
frame field results, but only hex meshes that was the main focus 


of their work. Therefore, we implemented both methods; there are 
few points worth noting: 

Sampling: In previous works the frame fields were sampled ei¬ 
ther on each tet face or on each tet. Instead we sample it on vertices, 
otherwise we would not be able to compare corresponding energies. 

Gimbal Lock: Both Huang and Li use Euler angles as variables 
in their L-BFGS optimization, which have numerical issues when 
the angles are close to the gimbal lock. Note that each frame can 
be represented by 48 triplets of equivalent Euler angles. In our im¬ 
plementation we choose the triplet maximizing the distance to the 
nearest gimbal lock. 

Rendering: Eor rendering purposes, we rely on a combination of 
techniques (Figure to show the spherical harmonics field, the 
frame field (locally and globally) and the field topology. 

3.6.1 Comparison with Huang's method. Recall that Huang et 
al proposed a method in two steps: 

—find an initial frame field by solving a linear system and project¬ 
ing the solution onto the manifold of feasible solutions 

—represent each frame by a triplet of Euler angles and optimize 

the smoothness using an L-BFGS descent method. 

Our impl ementation produc es results very similar to those pre¬ 
sented in [Huang et al. 20II| , but with significantly better timings. 
For example, the rockarm (Figure [TTJ with one million tetrahedra 
takes about 10 minutes on a single thread application on a Dell 
M6600 laptop compared to 155 minutes obtained by Huang et al 
on a two-thread i7 processor. 

Huang’s initialization is very similar to ours, however their 
boundary condition is not sufficient in this case (it requires the SH 
coefficients to be feasible). Moreover, they enforce the boundary 
condition with a penalty term that is very light (10“^ weight). As a 
consequence, their initial fields are almost constant everywhere (it 
maximizes the smoothness), with a topology very far from being 
optimal. The smoothing iterations are performed with much higher 
weight (10^) of the penalty term using L-BFGS. 

After the initialization step we measured the deviation of the 
field from the given constraints on the rockarm model. Note that the 
penalty term being the sum of deviations over all vertices, we can 
conclude that deviation at a given vertex belongs to [0, 2^7/12]. 
Thus on the rockarm the initial frame field has the average devia¬ 
tion of 0.34, whereas the maximum frame deviation is 0.96. 

If we use a much higher penalty weight to enforce the boundary 
constraint, we obtain an initial frame field with average deviation 
from constraints equal to 0.07 and maximum deviation equal to 
0.75. The field has better topology and L-BFGS converges faster 
for this initialization. The initialization provided by our method has 
average deviation from constraints equal to 10“® with maximum 
deviation of 10“^. 

Figure m gives an illustration, it compares three methods: 
Huang’s algorithm (left image) Huang’s algorithm with much 
higher penalty weight (middle image) and our initialization fol¬ 
lowed by Huang’s smoothing iterations (right). 

Huang’s paper is the pioneer work and the main contribution 
in [Huang et al. 2011) was the introduction of the energy used 
in frame field optimization, however the initialization is not very 
good and s moothing stage was also outperformed by later works 
(see §3.6.3| ). 

3.6.2 Comparison with Li’s method. Li’s initialization com¬ 
putes a 2D frame field on the surface, then propagates it inside 
the volume by advancing front. As a consequence, the initial field 
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Fig. 11. Comparison of initialization of Huang’s algorithm with their 
weight of boundary penalty term (10“^) (left/blue), with a much higher 
weight (10^) (middle/yellow), and our initialization (right/red). All 
smoothing iterations are performed by Huang’s algorithm. For the rocker 
model, the energies (we take the run with our initialization for 100%) are 
respectively 94, 7%, 101, 6% and (obviously) 100%. For the tet model, we 
obtain 91%, 89, 7% and 100%. 


perfectly matches boundary constraints, but is discontinuous across 
its medial axis. 

The smoothing iterations are performed by L-BFGS. It acts on a 
new set of variables that characterizes, per vertex, the rotation that 
brings the reference frame to the current frame. For the frames lo¬ 
cated inside the object, variables are the Euler angles as in Huang’s 
method. For frames located on the object boundary, the rotation is 
characterized by a single rotation angle around the normal vector. 

The frame field results presented in their paper were limited to 
an ellipsoid and a sphere (frame field design was not the primary 
objective). We guess that most of presented results were not fully 
automatically generated frame fields, as they wrote: “For instance, 
we use guiding boxes to modify the frames inside the narrow ears 
of Bunny (Figure 13-a) and the head of rock arm (Figure 13-b) to 
reduce singularities”. 

Moreover, before implementing the method, we thought that 
their algorithm was strongly limited by the original field topology 
from the sentence: “However, our propagation-based frame field 
initialization likely generates singular edges around the medial axis 
of the volume, and most of them cannot be eliminated by frame opti¬ 
mization.” Surprisingly, our implementation of their method is able 
to generate smooth frame fields automatically in most situations, 
even when the topology of the initial field is complex close to the 
medial axis. Our implementation is slightly different: 



Fig. 12. Li’s algorithm (red) is compared to our algorithm (green) on a 
one-finger bowling ball. The hole has a huge impact on the initialization due 
to the advancing front approach (second column, inside the yellow box). 
As a result, their initialization provides a field with a poor topology and 
smoothing iterations are not sufficient to find the expected topology (like 
ours). Our final energy is 86% of theirs. 



Fig. 13. Comparison of two fields generated by Li et al smoothing itera¬ 
tions: using their initialization (red) or ours (blue). Our energy in percent of 
theirs is 99%, 99.87%, 100.5%, 99.97%, 99%, 99.6%. The difference is 
always very low (< 1%) but always in our advantage except for the third 
model. 


—we initialize the 2D field by our 3D algorithm restricted to 
boundary vertices 

—we sample the field on vertices 

—and we prevent gimbal locks by a proper initialization of Euler 
angles. 

The only real failure case we found using their method is due to 
the front propagation algorithm: when a boundary frame is copied 
to a large number on inner samples. In this case, the L-BEGS solver 
is sometimes locked with a bad topology (see Eig.p^. 

On more complex examples, we have compared their algorithm 
against our initialization followed by their smoothing iterations. In 
most cases, we obtain an energy that is a bit better (Eig[^. We 
also observed that their topology often differs from ours (Eigp^, 
so we conjecture that our topology is somehow better. However, 
the quality of the field topology depends on the application, and is 
not well evaluated by the energy, even for topologies very far from 
being optimal (see e.g. Fig[n). 

3.6.3 Comparison of smoothing iterations. In the previous 
sections we have shown (Eig. m and that our method pro- 
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Fig. 14. Comparison of two fields generated by Li’s smoothing iterations: 
using their initialization (red) or ours (blue). Our energy in percent of theirs 
is 98.5%. Close-ups show the field where the singularity graphs diverge: 
one is inside the volume (leftmost) and others are on the object boundary. 
Our results are on the top row and theirs are on the bottom row with singu¬ 
larity encircled in white. 

vides the best initialization, however our smoothing iterations are 
not clearly better than others. 

Figure shows a comparison of three different smoothing 
strategies (our linearization, L-BFGS proposed by Huang et al and 
L-BFGS proposed by Li et al). In all three cases we use our method 
to initialize the field. L-BFGS smoothing proposed by Huang et al 
is the slowest in all test cases. First of all, in our implementation the 
time to evaluate the energy and the gradient is four times slower 
for the method by Huang et al than for the method by Li et al 
Second, the cruicial difference between these two methods is the 
way to enforce the boundary alignment: Huang et al use a penalty 
term, whereas Li et al use the set of variables directly satisfying the 
boundary constraints. In our test we noticed that usage of penalty 
terms increases the number of iterations to converge and interferes 
with topological choices to be made, leading to inferior final fields. 

We also noticed that the behaviour of linearization changes in 
function of how far the initialization is from the final solution. 

—First row of figure shows a simple case without topology 
changes, the smoothing iterations change the field geometry 
only. In this case the linearization of feasibility constraints works 
flawlessly, in two iterations the method converges, taking about 
the same time as the smoothing by Li et al 

—Middle row shows a second case, where few topology changes 
must be made. It slowes the linearization down, even if two iter¬ 
ations produce a reasonably good field. 

—Finally, the bottom row shows the case where the initial topol¬ 
ogy is really bad. The linearization method fails on this model: 
two first iterations are still very far from the final solution and 
to reach the minimum it requires four times more time than the 
method by Li et al 

We conclude that the best solution is our initialization followed 
by the optimization of Li et al. In practice, the implementation 
of our linearization smoothing iterations is almost free (incremen¬ 
tal with respect to our initialization algorithm), whereas Li et al 
smoothing algorithm is more difficult to implement. Moreover, in 
most cases few iterations of linearization steps suffice to obtain a 
fairly good field. As a consequence, we suggest starting with our 
smoothing iterations (almost free to implement), then possibly re¬ 
place it with Li et al smoothing algorithm if performances are not 
sufficient. 



Fig. 15. Comparison of smoothing iteration algorithms combined with our 
initialization algorithm. We compare Huang’s method (blue), Li’s method 
(green), our method (orange), and our method limited to two iterations (red). 
Top row compares the initialization (left) with the field after two iterations 
of our algorithm (right). Middle and bottom rows compare the field with two 
iterations of our smoothing algorithm with other smoothing strategies. Sin¬ 
gularity graphs reflect nicely the convergence of thee algorithms. We obtain 
energy (we take Li’s result for 100%) of resp. 99.98%, 100%, 99.98% 
and 100.5% on the sector, 101.5%, 100%, 100.4%, and 106.5% on the 
fandisk, and 116%, 100%, 109%, 185% on the one-finger bowling ball. 


Future works 

From an application point of view, our method makes it fast and 
easy to produce smooth frame fields. The smoothness is not neces¬ 
sary the optimal objective for applications such as hex remeshing 
(see Fig. [^, but can be used as a regulation term for more com¬ 
plex energies. It is also very easy to modify our method to add con¬ 
straints inside a volume. Figureshows a frame field constrained 
by faults in geological data. Note that the field is not constrained 
by the boundary of the model. Such a frame field is useful to steer 
element placement for fluid simulation used in oil exploration. 

From a theoritical point of view, it would be interesting to better 
understand the shape of the feasible set of a* (the 3D manifold 
embedded in 9L>). It could help to find a better projection algorithm 
than our current gradient descent. 
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RU^/2) = 


0 0 0 0 0 VT4/4 0 

0 -3/4 0 Vl/A 0 0 0 

0 0 0 0 0 72/4 0 

0 77/4 0 3/4 0 0 0 

0 0 0 0 3/8 0 

-714/4 0 -72/4 0 0 0 

0 0 0 0 75/4 0 

72/4 0 -714/4 0 0 0 

0 0 0 0 7^/8 0 -77/4 0 


-72/4 0 


714/4 0 

0 0 
0 7^/8 


0 

-77/4 

0 

1/8 . 


RW = R%{^m X RhiP) X R%{n/2y 


Fig. 20. The thin object (left) frame field is suitable to produce a hex mesh: 
its singularity graph is basically two singularities of index 1/4 and —1/4 
extruded in the 2 ; direction. The frame field of fat object (right) has a singu¬ 
larity graph that does not correspond to a hex mesh. 

Conclusion 

This work unifies the frame field design problem in 2D and 3D. 
Both problems are formulated with a similar representation of 
frames, constraints and objective function. As a consequence, they 
can also be solved by similar algorithms. 

From this analysis, we discovered that the best actual solution 
to produce smooth 3D frame fields is to initialize it by our pro¬ 
posed extens ion of [Kowalski et al. 2012| , followed by smoothing 
iterations of |Li et al. 2012| . The main drawback of this solution is 
requirement to implement two very different approaches (a sparse 
linear system solver a nd a L-BFGS descent) . A fair alternative is to 
use our extension of | |Kowalski et al. 20T^ , it is simple to imple¬ 
ment and requires a linear system solver only. In practice for our 
models we perform only two or three linearization iterations, how¬ 
ever if the initialization is a bad guess (e.g. the sphere), it can be 
insufficient. With this approach we are able to generate (on a lap¬ 
top) fields on the models up to few millions tetrahedra in less than 
10 minutes, refer to Figure[^for an illustration. 


RUa) = Rli^/^V X RUa) X i?^(7r/2) 

These matrices are calle d Wigner D-matrices and the literature 
on their construction is vast |Collado et al. 1989l|Blanco et al. 1997| 

|Choi et al. 1999[|Ivanic and Ruedenberg 1996| . However, as we are 

using degree 4 harmonics only, we performed the symbolic compu¬ 
tation as follows. 

The matrix R%{'y) is easy to compute, given a spherical har¬ 
monic F4,m(^5 0) we can rotate it around the z-axis by the angle 7 
by evaluating l 4 ,m 0 + 7 )- Then the element {i,j) of the matrix 

TT 27r 

is simply J f Y 4 ^i_ 5 ( 0 , (j)) • 0 + 7) sinOdOdcj). 

0 0 

Matrices R^^ (/3) and R% (a) are trickier to get, however we can 
use the fact that arbitrary rotation around the y-axis can be decom¬ 
posed into the rotation around x-axis by 7r/2 followed by a rotation 
around z-axis. To rotate a spherical harmonics l 4 ,m (^5 0) around 
the x-axis by 7r/2 we can perform the following substitution: 

{0,(1)) ^ (arccos(—sin(^) sin(0)), atan2(cos(^), sin(^) cos(0))). 

Then the matrix R%{7r/2) is calculated by evaluating double inte¬ 
grals of all possible products between basis functions and rotated 
functions. 


APPENDIX 


A.2 Linearization 


In order to linearize the constraints in our linear solver we define 
A. SH COOKBOOK matrices and E^ as follows: 

A.1 9D rotation 


Let us denote by R^, R^ and R^ 3 x 3 matrices of rotation around 
axes X, y and z respectively. Any frame / can be obtained as a 
composed rotation of the reference frame /, where the reference 
frame is the set of 6 unit vectors aligned with coordinate axes: 

f = R^{a)xRy{l3)xR^{j)xf. 




roooooo 0 -72 0 ■ 

00 000 0 -7772 0_-72 

00000 -3/72 0 -Yfn 0 

0 0 -TTO 0 -2,/V2 0 0 


0 0 
0 0 


0 TTo 0 

0 0_3/72 0 0 

0 77/2 0 3/72 0 
72 0 77/2 0 0 

L 0 72 0 0 0 


If {a, j3, 7 ) are Euler angles of rotation between a frame / and /, 
the representation vector a is calculated as a = R% (a) x R^ {{3) x 
R%{'y) X a, where R%, R^ and are 9 x 9 matrices of rotation 
defined as follows: 


^b{i) — 


El = 


- 0 72 0_ 

-72 0 7772 

0 -77^ 0 

0 0 -3/72 

0 0 0 
0 0 0 
0 0 0 
0 0 0 

- 0 0 0 


0 0 0 0 

0 0 0 0 

3/72 00 0 

0 0 0 0 

0 0 -TTo 0 

0 TTo 0 -3/72 

0 0 3/72 0_ 

0 0 0 77/2 

0 0 0 0 


0 0 - 
0 0 
0 0 
0 0 
0 0 
0 0 
-7772 0 
0 -72 

72 0 - 


” cos(47) 

0 

0 

0 

0 

0 

0 

0 

sin(47)" 










0 

cos(37) 

0 

0 

0 

0 

0 

sin(37) 

0 


0 

0 

0 

0 

0 

0 

0 

0 

4" 

0 

0 

cos( 27 ) 

0 

0 

0 

sin( 27 ) 

0 

0 


0 

0 

0 

0 

0 

0 

0 

3 

0 

0 

0 

0 

cos('y) Osinfo') 

0 

0 

0 


0 

0 

0 

0 

0 

0 

2 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 


0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

- sin( 

7 

0 

0 

7 

0 

0 

0 

Eb = 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

- sin( 27 ) 

0 

0 

0 

0 

0 

to 

0 

0 


0 

0 

0 

-1 

0 

0 

0 

0 

0 

0 

- sin(37) 

0 

0 

0 

0 

0 

cos(37) 

0 


0 

0 

-2 

0 

0 

0 

0 

0 

0 

— sin(47) 

0 

0 

0 

0 

0 

0 

0 

cos(47) 


0 

-3 

0 

0 

0 

0 

0 

0 

0 










-4 

0 

0 

0 

0 

0 

0 

0 

0 
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It is easy to verify that these matrices are chosen to verify the 

following equations for small rotations a, 7 : 

RUo^) = ^ aE% ^ o{\a\) 

Rim = 19.9+I 3 EI+om) 

Rhm = i9.9+iEi+o{\m- 

Finally, for small rotations the multiplication is commutative: 

RBia,0,j) = Rim X Rim X Rim = 

— ^9x9 + + l/^l + ItD- 
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Fig. 16. Our results are shown using combinations of the following ren¬ 
dering techniques. We can plot for each vertex i its F'^ (upper left), or its 
frame as a cube(upper right). We can show the singular lets (lower left), 
or a smoothed and refined version (lower middle) to better see it in 3D 
(thanks to the lighting). The field inside the volume can be rendered by 
curved french fries. 



Fig. 17. The initialization of |Huang et al. 2011] (left) is a constant 
frame field whereas ours (right) is aligned to the boundary. Their F* are 
all equal, and very far from being feasible, making it possible to violate 
the boundary condition. 



Fig. 18. Results on CAD objects. Names are (from left to right): Anc, Crank, 40head. 



Fig. 19. 3D frame field constrained by faults in geological data. 
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